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A finite element formulation with linear triangular elements was used to 
solve the steady-state, two-dimensional conduction heat transfer equation 
in the condenser wall section of an internally finned rotating heat pipe. 

A FORTRAN program using this method was coupled with the ADS program 

for automated design of the internal heat pipe fin geometry to optimize 
heat transfer. An increase in surface area, which increases heat transfer, 
also increases the condensate level, which decreases heat transfer. 

The additional condensate level does not offset the advantage gained by 
the increased surface area. The investigation provided combinations of 
fin half angle, number of fins, and fin height for an optimum design. 


Water is used as the working fluid and the heat pipe is constructed from 
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I. INTRODUCTION 
A. THE ROTATING HEAT PIPE 

The rotating heat pipe is a closed container designed to transfer a 
large amount of heat in rotating machinery. Since the heat pipe operates 
on a closed two-phase cycle, the heat transfer capacity is greater than for 
solid conductors. Also, the thermal response time is less than with solid 
conductors. The three major elemental parts of the rotating heat pipe are: 
a cylindrical evaporator, a truncated cone condenser, and a fixed amount 
of working fluid as shown in figure l. 

An annulus is formed by the working fluid in the evaporator. This 
occurs at rotationary speeds above the critical speed. The addition of heat 
to the evaporator vaporizes the working fluid. A pressure differential 
between the evaporator and the condenser causes the vapor to flow 
towards the condenser. The vapor is transported to the condenser with its 
latent heat of vaporization. Condensation of the vapor on the inner wall is 
caused by external cooling. This condensation releases the latent heat of 
evaporation. This condensate is forced to flow back to the evaporator by 
a component, acting along the condenser wall, of the centrifugal force 
which is caused by the rotation of the heat pipe. As the condensate 
collects in the evaporator the cycle is repeated. 

Since the evaporator and condenser portions of a heat pipe function 
independently, needing only common liquid and vapor streams, the area 


over which heat is introduced can differ in size and shape from the area 


heat out 





Fizure Ll. Schematic Drawing of 4 Rotating deat Pipe 


Z 


N 


over which it is rejected, provided that the rate at which the lquid is 
vaporized does not exceed the rate at which it can be condensed. 
Therefore, high heat fluxes generated over relatively small areas can be 
dissipated over larger areas with reduced heat fluxes; allowing a 
cylindrical evaporator and a truncated cone condenser. 

Capillary action acts to drive the condensate back to the evaporator 
in a conventional heat pipe. No limitation due to capillary action is 
encountered in a rotating heat pipe nor are external pumps or gravity 
depended on for the flow of the working fluid. Therefore, the rotating heat 


pipe can be used in any orientation [Ref. 1]. 


B. OPERATING LIMITS OF A ROTATING HEAT PIPE 

The first theoretical investigation of the rotating heat pipe conducted 
at the Naval Postgraduate School was performed by Ballback [Ref. 2] in 
1969. Various fluid dynamic mechanisms limit the performance of a rotating 
heat pipe. Ballback [Ref. 2] studied these mechanisms and an estimation of 
the sonic limit, boiling limit, entrainment limit, and condensing limit of 
performance was made. 

1. The Sonic Limit 

The maximum flow of the — is set by the choked flow 

condition in the rotating heat pipe. This limiting vapor flow rate occurs 
when the heat flux is increased and limits the amount of energy the vapor 
can eee The rotating heat pipe effectiveness is reduced due to this 


limitation. The limiting heat transfer rate due to this condition is: 


Q, = p,U Aly m 


and the vapor velocity is considered to be sonic, 


(2) 


tN | 


U, = c = (¢,kRT) 


where 


U, = velocity of the vapor in ft/sec, and 


A = cross sectional area for the vapor flow in ft 


sonic velocity in ft/sec 


S 
Gy = gravitational constant 


k = ratio of specific heats 


R = gas constant in ft-lbf/lbm R 


ay 


absolute temperature in degrees Rankine 
py = density of vapor in lbm/ft° 
2. The Boiling Limit 
The transition from nucleate to film boiling was hypothesized by 
Kutateladze [Ref. 3] to be a completely hydrodynamic process. He 
determined the following theoretical formula for predicting the burnout 


Elise 


p|—m 


( 
Q, - K/p, A,h,log(p -p,)} 7 


where 
K = constant value 


A, = heat transfer area in the boiler in et? 


heg = latent heat vaporization in Btu/lbm 


Oo = surface tension in lbf/ft 

g = acceleration of gravity in ft/ hr? 
oc = density of fluid in lbm/ft° 

py = density of vapor in lbm/ft?. 


A constant value for K in the range of 0.13 to 0.19 was suggested by the 
experimental data obtained by Kutateladze. 
3. The Entrainment Limit 
The flooding constraint in a wickless heat pipe was examined by 


Sakhuja [Ref. 4]. The correlation he developed is: 


A,C*h,./gD(p -p)P, 
Q 


at 
Li+(p Jp y*P 


(4) 


where 


tO 
be 
\ 


= heat transfer rate in Btu/hr 


= flow rate in ft? 


Hw - 
pa 
! 


C = dimensionless constant, 0.725 for tube with sharp edged flange 


he, = latent heat of vaporization in Btu/lbm 


g = acceleration due to gravity in ft/hr? 
D = inside diameter of heat oe Kel ome 

o¢ = density of the fluid in lbm/ft° 

py = density of the vapor in lbm/ft?. 


4. The Condensing Limit 
The condenser section of a rotating heat pipe was modeled as a 
truncated cone by Ballhack [Ref. 2]. Using this model, the condensation 
limitation for a rotating heat pipe was determined by Ballback [Ref. 2]. He 


developed the following condensation limit: 


] 
2 Kp a*hy(T,-T.)” i 


[(R,+Lsino)*?-Ro°}* () 
3 in? ° 
p sin” 


Q, = nt 


where 
Qt = total heat transfer rate in Btu/hr 
ks = thermal conductivity of the condensate film in Btu/hr-ft-F 
os = density of fluid in lbm/ft° 


@ = angular velocity in l/hr 

he, = latent heat of vaporization in Btu/lbm 
T; = saturation temperature in F 

T, = inside wall temperature in F 


viscosity of fluid in lbm/ft-hr 


: = 
h 
it 


@ = half cone angle in degrees 
Ry = minimum wall radius in ft . 
L = length along the wall of the condenser in feet. 
The geometry and speed of the rotating heat pipe, and the physical 
properties of the working fluid comprise the condensing lmit equation. 
For a rotating heat pipe with the physical characteristics as 
shown in Table I, the amount of heat that can be transferred from the 


rotating heat pipe is limited by the condensing limit. However, the 


limitations imposed by the sonic limit, boiling limit, and entrainment limit 
may become important as the heat pipe geometry and operating conditions 
vary. 


TABLE I. ROTATING HEAT PIPE SPECIFICATIONS 
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| Length 9.000 inches | 


f 


| Minimum Diameter 1.55 inches 
| Wall Thickness 0.03125 inches | 


Internal Half Angle 1.000 degree 
| Rotating Speed 3600 RPM 


To enhance the heat transfer capacity of the rotating heat pipe, 


ey Se 


internally finned condensers have been used to raise the condensing limit 
line. Thinner films occur near the ridges of the fins while thicker films 
occur in the troughs. The thinner film on the ridges provides a lower 
thermal resistance to heat flow, while the thicker film in the trough 
provides a higher resistance. A compromise between the improvement. on 
the ridges and the degradation in the troughs is necessary for an overall 


heat transfer improvement [Ref. 5]. 


C. ANALYSIS OF THE INTERNALLY FINNED ROTATING HEAT PIPE | 
Schafer [Ref. 6] developed an analytical model for a meat pipe ae 
a triangular fin profile (figure 2). This model was developed in order toa 
raise the condensing limit by the addition of internal fins. An assumption 
of one-dimensional heat conduction through the wall and fin was made for 
Schafer’s model. 
A two-dimensional heat conduction model using a Finite Element Method 


was developed for this same case by Corley [Ref. 7]. A parabolc 


temperature distribution along the fin surface was assumed hy Corley [Ref. 
7]. A significant improvement in the predicted heat transfer performance 
was indicated by his results. By using the two-dimensional model an 
increase of approximately 75 percent in the heat transfer performance was 
seen over the results from the use of the one dimensional model. However, 
Corley (Ref. 7] noted that at the fin apex an error as great as 50 percent 
was possible which could result in the total heat transfer being in error 
as much as 15 percent. 

A modification was made to Corley's computer program by Tantrakul 
[Ref. 8]. In order to minimize the heat transfer error at the apex of the 
fin Tantrakul increased the number of finite elements used. His results 
with this modification converged with the results of Corley. 

Purnomo developed a linear triangular finite element model (figure 3) 
used in a two-dimensional. Finite Element Method solution. Purnomo's [Ref. 
1] Finite Element Method program also worked and converged. To maximize 
the heat transfer from the rotating heat pipe the condenser geometry was 
varied. Using Purnomo's code parametric studies were conducted. However, 
the best geometry was not indicated in these studies. Purnomo's code was 
written to perform one analysis at a time. Davis [Ref. 9] modified 
Purnomo’s code to allow for numerous analysis to be made using the 
optimization code COPES/CONMIN. Davis’ Finite Element .Method code 
incorporating the optimization worked and converged, resulting in an 


optimum design for an internally finned rotating heat pipe. 





Figure 2. Internally Finned Condenser Geometry, Showing Fins, 
Troughs and Lines of Symmetry. 





tc @ thickness 


b® height of the fin 


$ = condenser nalf angle 
a half of the fia width 
c= half of the srough width 


Figure 3. Condenser Geometry Considered with 40 Linear Triangular 
Finite Elements. 


D. THESIS OBJECTIVES 
The objectives of this thesis were: 
1. To modify Davis’ [Ref. 9] computer program so that it is compatible 
with the ADS (Automated Design Synthesis) program [Ref. 10] and 


can be used for analysis and automated design of rotating heat 
plpes., 
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To use the resulting program to obtain an optimum design for an 
internally finned rotating heat pipe to obtain experimental data to 
compare with the analytical results. 


To use the resulting program to obtain numerical results in place of 
data obtained from costly experimental operations. 


tal 


II. NUMERICAL OPTIMIZATION 
A. BACKGROUND 

The parameter that is minimized or maximized during the design 
process is called the design objective. The design objective is minimized 
or maximized by changing the design variables within the design constraint 
limitations. This process is called numerical optimization. An assortment of 
physical, aesthetic, economic and, on occasion, political limitations must he 
met by the design constraints for the design to be acceptable. For the 
optimization process to work, the design criteria must be described in 
numerical terms. This is not always easy. 

A computer program a be written to perform tedious and repetitive 
calculations necessary to optimize the problem once it is stated in 
numerical terms. For this reason, computer analysis 1s commonplace in most 
engineering organizations. For example, in heat transfer design the 
configuration, materials, and method of heat removal may he defined and 
a finite element analysis computer code is used to calculate temperatures, 
heat transfer rates, and other response quantities of interest. If any of 
these parameters are not within prescribed bounds, the engineer may 
change the method of cooling or other defined quantity and rerun the 
program. The engineer makes the actual design decisions, the computer 
code only provides the analysis of a proposed design. This is the commonly 


used approach which is called computer-aided design. 
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Analysis codes are commonly used for tradeoff studies. For example, 
an analysis code might be run on the distance a truck can go on a tank 
of fuel. For different loads, different distances are calculated which can 
be used in a range-payload study. 

Fully automated design is the logical next step to computer-aided 
design. The computer makes the actual design decisions or trade-off 
studies based on input criteria in fully automated design. Minimal 
information is requested from the operator during the actual design 
process. Numerical optimization offers numerous improvements over the 
traditional approach to design. These improvements include: time reduction 
in design decision making; a rational, directed design procedure; and the 
procedure is unbiased by intuition or experience. The probability of 
obtaining a non-traditional solution is thereby improved. Engineering 
intuition and experience are still necessary to decide if the design obtained 


is an improvement and feasible. 


B. AUTOMATED DESIGN SYNTHESIS (ADS) 

Vanderplaats [Ref. 10] developed a general purpose numerical 
optimization program containing a variety of algorithms, ADS. ADS is a 
FORTRAN program that optimizes a numerically defined objective function 
subject to a set of constraint limits. The solution of the problem is 
separated into three levels: 

1. Strategy - Optimization strategy such as Augmented Lagrange 
Multiplier method or Sequential Linear Programming. 
Ze ; Optimizer - Actual algorithm to perform the optimization. 


3. One-Dimensional Search - Line search routine used hy optimizer. 


Jes 


Flexibility to solve a wide variety of engineering design problems is 
given by the combinations of nine strategies, five optimizers, and eight 
one-dimensional search options. The following definitions are necessary to 
discuss the use of ADS: 

1. Design Variables - Those parameters which the optimization program 
is permitted to change within allowed bounds in order to improve 
the design. Design variables appear only on the right hand side of 
an equation and are continuous. 

2. Design Constraints - An inequality constraint requires that some 
function of the design variable(s) remain less than a specified value. 
Design constraints may be linear or nonlinear, implicit or explicit, 
but they must be continuous functions of the design variable. 

3. Objective Function - The parameter which is going to be minimized 
or maximized during the optimization process. The objective function 
may be linear or nonlinear, implicit or explicit, and must be a 
continuous function of the design variables. The objective function 
usually appear on the left side of an equation. 

Cz PROGRAMMING GUIDELINES 

Any computer code developed for engineering analysis should be 
written in such a way that it is easily coupled to a general purpose 
optimization program such as ADS. Therefore, a general programming 
practice is outlined here which in no way inhibits the use of the computer 
program in its traditional role as an analytical tool, but allows for simple 
adaption to ADS. 

ADS is called by a user-supplied calling program. ADS does not. call 
any user-supplied subroutines. Instead, ADS returns control to the calling 
program when function or gradient information is needed. The required 


information is evaluated and ADS is called again. This provides considerable 


flexibility in program organization and restart capabilities. Various internal 
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parameters are defined on the first call to ADS which work well for the 
"average'’’ optimization task. However, it is often desirable to change these 
in order to gain maximum utility of the program. Figure 4 is the program 
flow diagram for the case where the user wishes to over-ride one or more 
internal parameters, such as scaling, convergence criteria, or maximum 
number of iterations. 

After initialization of basic parameters and arrays, the information 
parameter, INFO, is set to -2. ADS is then called to initialize all internal 
Parameters and allocate storage space for internal arrays. Control is then 
returned to the user, at which point these parameters, for example 
convergence criteria, can be overridden if desired. At this point, INFO will 
have a value of 1 and the user must evaluate the objective function, OBJ, 
and constraint functions. ADS is called again and the optimization proceeds. 
Since, in this case, the gradient calculation control, IGRAD, has a value of 
zero, all gradient information is calculated by finite difference within ADS. 


When INFO has a value of zero, optimization is complete. 


Us 


BEGIN 

DIMENSION ARRAYS 

DEFINE BASIC VARIABLES 
IGRAD=0 (USE FINITE DIFFERENCE GRADIENTS) 
INFO = ~2 
CALL ADS (INFO...) 
IF INFO = 0, EXIT. ERROR WAS DETECTED 

ELSE 


OVER-RIDE DEFAULT PARAMETERS IN ARRAYS WK AND IWK IF 
DESIRED 


CALL ADS (INFO...) 


NO YES 
INFO = 0 
EVALUATE OBJECTIVE EXIT OPTIMIZATION 
AND CONSTRAINTS IS COMPLETE 
Figure 4. Program Flow Logic: Over-Ride Default Parameters, Finite 


Difference Gradients [Ref. 10]. 


16 


A. 


Ill. FINITE ELEMENT SOLUTION 


REVIEW OF THE PREVIOUS ANALYSIS 


As stated previously, the heat transfer solution for a one-dimensional 


model of an internally finned rotating heat pipe was studied by Schafer 


[Ref. 6]. The two-dimensional model was studied by Corley [Ref. 7]. The 


same assumptions and boundary conditions, similar to those used in the 


Nusslet analysis of film condensation on a flat wall, and based upon the 


analysis of Ballback [Ref. 2] were used for both. The more important of 


these assumptions are: 


Sie 


steady state operation, 
film condensation, as opposed to dropwise condensation, 


laminar flow of the condensate film along both the fin and the 
trough, 


static balance of forces within the condensate, 


one-dimensional conduction heat transfer through the falm thickness 
(no convective heat transfer in the condensate film), 


no liquid-vapor interfacial shear forces, 
no condensate subcooling, 


zero heat flux boundary conditions on both sides of the wall section 
(symmetry conditions), as shown in figure 5, 


saturation temperature at the fin apex, 


10. zero film thickness at the fin apex, and 


ll. negligible curvature of the condenser wall. 


i 


Figure 3 shows the linear triangular finite element model developed 
by Purnomo [Ref. 1] for use in obtaining a two-dimensional Finite Element 
solution. 

The assumption that was used by Corley [Ref. 7] that the fin apex was 
at the saturation temperature of the working fluid was modified by 
Purnomo [Ref. 1]. The value of the temperature at the apex of the fin was 
allowed to float and a parabolic temperature distribution was assumed along 
the fin surface. 

Purnomo's problem statement for the formulation of the Finite Element 


Method as shown in figure 5 is: 
es Se (6) 
with the following -boundary conditions: 
1. along boundary 1, -K o7/on-h,(T-Tsat) 
2 meee boundary 3, -K oT/on=h,(T-T~) 


oF. 


On 


3. along boundaries 2 and 4, 


A detailed description of the numerical formulation is presented in his 
thesis. 
Davis used Constrained Function Minimization (CONMIN) as an 


optimization program. CONMIN is a FORTRAN program in subroutine form. 
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Differential Equation and Boundary Conditions Considered 


ysis of Purnomo [Ref. 1]. 


Me, 


Vanderplaats [Ref. 11] developed the Control Program for Engineering 
Synthesis (COPES) as a main program to simplify the use of CONMIN. Davis’ 
computer program was written in subroutine form with SUBROUTINE ANALIZ 
(ICALC) as the main routine. The name ANALIZ is compatible with the COPES 
program and ICALC is a calculation control. Subroutine ANALIZ calls other 
subroutines as needed: 

Qe the pepe "COORD" used to define positions of system coordinate 

poln 


2. the routine "FORMAF” used to formulate the Finite Element Method 
equations, 


3. the routine "BANDEC”™ as an equation solver for a symmetric matrix 
which has been transformed into banded form, and 


4. the routine "DPLORT” used to compute the roots of a real polynomial 
using a Newton-Raphson derivative technique. 

B. THE COMPUTER PROGRAM DEVELOPMENT 
The basis for the present analysis code is Davis’ [Ref. 9] two- 
dimensional finite element program. Davis’ code was checked for validity as 
the first task undertaken in the development of this thesis. An error was 
discovered in calculating the fin condensate film thickness (AZS). In the 
initial calculation of HDEN only, the cubed term was merely multiplied by 


three. The correct form of the equation is shown below: 


HDEN = -a,z3/3-b,z?/2+2(T,,-T,) 


The effect of this error was minimal since subsequent equations were 


correct, a 0.00016% difference in the condensate level was noted. 
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The next task undertaken was to adapt the analysis code to permit 
automated design analysis using ADS. Many modifications were made, some 
of which are mentioned here. The program was rewritten to include a main 
program from which ADS is called and subroutines to perform various 
mathematical functions. Subroutine ANALIZ was deleted since the double 
precision version of ADS (DADS) was used. Modifications were also made to 
generalize the code and minimize the changes needed when varying the 
number of finite elements used. 

A listing of the revised computer program is included as the 


Appendix. 


C. DESIGN OPTIMIZATION 

There are thirteen parameters that can be used as design variables. 
There are geometric or functional parameters of the rotating heat pipe or 
the properties of the working fluid or environment. The possible design 
variables, possible constraint functions, and the objective function are 
listed below in Fortran. This code can pursue a wide variety of design 
problems. 

The addition of fins by the designer increases the surface area which 
increases the heat transfer rate through the condenser wall. However, the 
addition of fins decreases the cross-sectional area through each fin for 
conduction and decreases the trough einen which increases the condensate 
thickness in the trough. The increased condensate thickness decreases the 
heat transfer and if increased to the point of covering the fins it could 


dramatically reduce the heat transfer through the fin. 


an 


TABLE II. DESIGN VARIABLES 


| CANGL (cone half angle) 
CLI (condenser length 

| FANGL (fin half angle) 

| HINF (external convective heat transfer coefficient) 

| 


| 
| BFIN (fin height) | 
| | 


| NFIN (number of fins) ! 
i R2I (intermediate radius) : 
i RBASEI (inside radius of condenser base) 
| RMP (rotational speed of the heat pipe) 
THICKI (condenser wall thickness) 
| TINE (external temperature) 
i TSS (saturation temperature) ) 
TZ (nodal point temperature) H 


| CONSTRAINT FUNCTIONS 


BOA (ratio of fin height over fin width) | 
ZOA (ratio of trough width over fin width) 
|; DEL(NT) (condensate thickness) | | 


| OBJECTIVE FUNCTION 
| OBJ 


== : = = = = = SS I | 


The purpose of the design study was to determine the fin height, 
number of fins, and fin half angle which would maximize the heat transfer 
rate. It was then decided that the design variables would be BFIN, FANGL, 
and NFIN. The number of fins was chosen vice the ratio of trough width 
to fin width (ZOA) as was previously used. This decision was made since 
it 1s easier to think in terms of the number of fins vice a ratio. Other 
potential design variables were held constant. The objective function to be 
maximized was OBJ=QT+QTF, the heat transfer through the fin plus the heat. 
transfer through the trough. 

The code was run with the three design variables using the internal 
scalar default parameters in ADS. The objective function was calculated 


using the input values of the design variables. The ADS program then 
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changed the first design variable keeping the other two design variables 
at the input values. The calculations were made aging yielding a different 
value for the objective function. The new value would then be compared 
to the previous value of the objective function. If the difference in the 
objective function was not greater then the internal scalar default value, 
the first design variable would be returned to its original value and the 
process repeated with the second design variable. If the difference in the 
objective function was still not areater then the internal scalar default 
value, the second design variable would be returned to its initial value and 
the process repeated with the third design variable. 

Each time the program was run, the optimization code would choose 
BFIN as the design variable to change first as it had the greater influence 
on the objective function. The remaining two variables would then either 
be kept constant or the number of fins would be maximized and the fin 
half angle minimized. Consistent results were not obtained with this method. 

To improve the results, the internal scalar parameters were modified. 
These modifications included the constraint tolerances, the absolute and 
relative convergence criteria, the absolute and relative change in the 
design variables, the absolute and relative change in the objective 
function, the minimum absolute value of the finite difference step when 
calculation gradients, and the initial relative move limit. Better results were 
obtained as seen in the higher value for the objective function. However, 
the results were still not consistent depending on the initial values used 
for the three design variables. At this point, it was decided to concentrate 
on one design variable and on the basis of the previous calculations the 


design variable chosen was BFIN. 


Ze 


The external surface temperature was set equal to the working fluid 
saturation temperature and the theoretical upper limit on the heat transfer 
was calculated for comparison. An assumption is made that there is no 
thermal resistance across the condensate or the condenser wall. The upper 
limit of the heat transfer rate was predicted to be 69,492 BTU/HR using the 


following formula: 


Qa ~ 2url(T,.7-T.) (7) 


where 


h = outside convective heat transfer coefficient (5000 BTU/HR‘FT:'F) 


average outside radius of condenser wall (0.07373 FT) 


rh 
ti 


1 = condenser length (0.75 FT) 

Tyai1 = temperature of the outside wall (100°F) 

T, = ambiont temperature (60°F) 
Based on engineering judgement certain constraints were placed on the 
design. These constraints were applied to the number of fins (not to 
exceed 400) and the minimum fin half angle (not to be less than 10 
degrees). The pales were based on structural and manufacturing 


considerations. 
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A. INTRODUCTION 
The purpose of the design optimization was to’ maximize the heat 
transfer rate. This was accomplished by using the computer code in 


conjunction with ADS. Numerical results are discussed below. 


B. CONSTRAINED OPTIMIZATION 

In the design problem undertaken to determine the optimum internal 
geometry for the maximum heat transfer, numerous runs were made for a 
condenser made of copper. This material has a thermal conductivity of 230 
BTU/HR'FT'F. The working fluid was water. : 

Since the fin height (BFIN) was the design variable, the initial runs 
investigated whether there was an optimum fin height. Initially the fin half 
angle was held constant at 10 degrees and the number of fins was varied 
from 150 to 400. In each case, for the optimum design, the fin height was 
maximized and the trough was eliminated, as seen in figures 6 and 7. 

The number of fins was then held constant and the fin half angle was 
varied from 10 to 25 degrees with the fin height remaining the design 
variable. Once again, the greatest heat transfer rate was achieved with the 
highest fin height for each number of fins (figure 8.) 

As seen in figure 9, the highest heat transfer rate achieved was for 


400 fins with a 10 degree fin half angle and a fin height. of 0.0345 inches. 
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Figure 6. Heat Transfer Rate vs. Fin Height. 


26 


(BTU/HR) 


HEAT TRANSFER RATE 


62500 
£ 6000C 
37500 


$500C 
52500 
$0000 
4750¢ 
4500C : 


42500 


40000 


37800 


L 
r 


Lo 
s=000 : ; 3 
0.005 0.010 90.015 9.020 90.025 0.030 0.035 0.040 0.045 0.050 


FIN HEIGHT (INCHES) 


LEGEND 
O = NUMBER OF FiINS=S00 
o = NUMEER OF FINS=SS0 
& = NUMEER OF FINS==CO 


Figure 7. Heat Transfer Rate vs. Fin Height. 


27 


HEAT TRANSFER RATE (BTU/HR) 
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Figure 8. Heat Transfer Rate vs. Fin Half Angle (Optimum Fin Height). 
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Figure 9. Heat Transfer Rate vs. Number of Fins (Optimum Design). 
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Figure 10 shows a plot of heat transfer rate versus fin half angle for 
a condenser with between 100 and 400 fins. The heat transfer rate, as a 
function of fin half angle for a constant fin height, increases as the fin 
half angle increases. Davis [Ref. 9] concluded that the heat transfer rate 
increased with an increase in fin half angle. This is correct if the fin 
height is kept constant. As the fin half angle increases, the surface area 
of the fin increases. The added surface area also has a thinner film of 
condensate on it which offers lower thermal resistance. The trough area 
decreases and the condensate film in the trough thickens, increasing the 
thermal resistance. This degradation does not offset the gain in the heat 
transfer rate caused by the fin. 

| For external heat transfer coefficients varying from 1000-50,000 
BTU/HReFT“ oF, the same optimum design geometry for a maximum heat 
transfer rate was obtained, which is stated in table III below. Figure 11 
shows the strong influence the external heat transfer coefficient has on 
the heat transfer rate. 

In figure 12, the effect of the rotating speed on the heat transfer 
rate is seen. As the rotational speed increases, the heat transfer rate 
increases, this 1s caused by an increase in the element heat transfer 
coefficient and a decrease in the condensate thickness on the fin which 


lowers the thermal resistance. 
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Heat Transfer Rate vs. Fin Half Angle (Constant Fin 
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HEAT TRANSFER RATE (BTU/HR) 
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Figure ll. Heat Transfer Rate vs. Heat Transfer Coefficient. 
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Figure 12. Heat Transfer Rate vs. Number of Fins (RPM Variation). 
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Figure 13. Condensate Level vs. Position (100-400 Fins). 
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eet = vx 





When the number of fins was increased from 100 to 400, maintaining 
the same fin half angle, the condensate level decreased with the increased 
heat transfer rate (figure 13). This occurred because the thinner film over 


the fins decreased the resistance across the film which raised the 


‘temperatures along the fin, which in combination with the lower height fins 


increased the temperature on the outside of the pipe. This increase in 
temperature brings the operation closer to the condensing limit. When the 
fin half angle was increased for a specified number and height of fins, the 
condensate level decreased due to the increased trough width (figure 14). 


TABLE III. OPTIMIZATION RESULTS FOR 
VARIOUS HEAT TRANSFER COEFFICIENTS 


& OPTIMUM DESIGN . 


Fin Height 0.0345 inches 
Fin Half Angle 10.0 degrees 
Number of Fins 400 


Nespertial HEAT TRANSFER RATE 


(BTU/HReFT“eP) (BTU/HR) 


13765 
62252 
261003 


Figures 15 and 16 show the effect the increase of surface area has 
on the heat transfer rate. In figure 15, the fin height for 400 fins with a 
10 degree half angle is plotted against the heat transfer rate. As the fin 
height is increased up to a maximum value of 0.0345 inches the resulting 
design is a sawtooth. Figure 16 shows the effect adding more fins has on 


the heat transfer rate for a constant height fin. The greatest increase in 
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the heat transfer rate is seen from the addition of 100 fins from a smooth 
tube. 

In figure 17, the ratio of the actual surface area over the surface 
area for a smooth tube is plotted versus the number of fins. The fin 
height and the fin half angle are both held constant. As expected, the 
surface area ratio increases in a linear manner as the number of fins 
increases. The ratio of the heat transfer rate over the heat transfer rate 
. for a smooth tube is plotted for two different heat trensfer coefficients, 
h=1000 BTU/HReFT“*F and 5000 BTU/HR*FT“eF. An increase in the number 
of fins results in not only an increase in the area ratio but also an 
increase in the heat transfer ratio. The increase in the heat transfer ratio 
is greatest when going from a smooth tube to a tube with 100 fins. The 
heat transfer ratio increase is greater for the heat transfer coefficient 
equal to 5000 BTU/HReFTePF. This is because the heat transfer coefficient 


has a direct effect on the heat transfer rate, that is, 
Q - hA(T-T ao) 


Both curves approach an asymptotic value. However, the curve with the 
lower heat transfer coefficient approaches this asymptotic value with a 
fewer number of fins. Additionally, in view of the relatively small increase 
in the heat transfer ratio by the addition of fins for the lower heat 
transfer coefficient, consideration should be given to the cost of 


manufacturing the fins versus the benefit derived by their addition. 
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CONDENSATE THICKNESS (FEET) 


Figure 14. Condensate Level vs. Position (400 Fins, 10-25 Degree Fin 


Half Angle). 
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Figure 15. Heat Transfer Rate vs. Fin Height (400 Fins). 
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Figure 17. Heat Transfer and Area Ratios vs. Number of Fins. 
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C. AUTOMATED DESIGN SYNTHESIS (ADS) 

This optimization project was done on the IBM mainframe using the 
ADS optimization program. As stated previously, ADS is a general purpose 
numerical optimization program with a variety of algorithms that can be 
used to tailor the solution. The solution is separated into three levels in 
ADS: strategy, optimizer, and one-dimensional search. For this problem the 


following combination of algorithms were used: 


1. Strategy: Sequential Linear Programming 


2. Optimizer: Modified Method of Feasible Directions for constrained 
minimization 


3. One-Dimensional Search: Golden Section Method followed by 
polynomial interpolation 


The strategy used linearizes a nonlinear problem by a first order 
Taylor series expansion of the objective and constraint functions. The 
solution to this linear approximation is obtained. The problem is linearized 
again about this point and the new problem is solved with the process 
being repeated until a precise solution is achieved. 

The optimizer chosen is used to find a search direction which will 
minimize the objective function while maintaining feasibility. 

The combination of strategy, optimizer, and one-dimensional search 
chosen is not the only one available, nor is it necessarily the most 
efficient. It did yield eee that were maximized and were within the 
constraint tolerances. 

The ADS optimization program is complicated by the numerous internal 


parameters which must be changed to obtain an optimal design. 
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Complications arise when there is a vast difference in the scales of the 
design variables. The design variables themselves must be scaled which is 
further complicated when the variable itself covers a wide range. ADS also, 
in certain cases, allows for constraints to be violated. In some instances 
this might be acceptable but not in this case. 

ADS does not have a scoping mechanism, that is the ability to 
decrease the rate of change of the design variable, and therefore to obtain 
a precise answer the internal parameters must be changed repeatedly. ADS 
also does not recognize integers, all numbers are real therefore, depending 


on the answer given, it may be necessary to round up or down. 
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V. CONCLUSIONS 


1. For an independent increase in fin half angle, rotational speed or 
number of fins an increase is seen in the heat transfer rate. As the 
parameters increase, the heat transfer rate levels off at the theoretical 
maximum heat transfer rate for the heat pipe. A decrease in the fin half 
angle with a mor restentinic increase in the fin height increases the heat 
transfer rate. If the fin half angle is decreased while the fin height is 
kept constant, then the heat transfer rate decreases. 

2. Maximum heat transfer occurs for the same fin geometry 
regardless of the external heat transfer Soastisiene. For a _ specific 
condenser radius, as many fins as possible should be machined with a 
minimum fin half angle at the maximum fin height. aes 

3. The computer code can be used for single analysis or the 
automated design of an internally finned rotating heat pipe. 

4. The benefit of adding fins is dependent on the external heat 
transfer coefficient. Consideration of the cost of manufacturing the fins 


versus the increase in the heat transfer rate should be made. 
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VI. RECOMMENDATIONS 


l. Analyze different shaped fins including rectangular and curved. 

2. Modify the code to allow for silmutaneous variations of more than 
one variable. 

3. Use different working fluids and heat pipe materials to see if a 
ene internal geometry occurs for the maximum heat transfer rate. 


4. Modify the code to use the DOT optimization program vice ADS. 
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APPENDIX: PROGRAM LISTING 


PROGRAM OLSON 
HK IRI I II II IIIT IIT I III TI ITT TTT TOT TOT TOIT TOI TOTTI IOI IOI OTIC Hk 


ke x 
* ANALYSIS OF ROTATING HEAT PIPE , USING TRIANGULAR “ 
* ELEMENT MODEL wd 
* COMPILED BY MAJOR IGNATIUS.S.PURNOMO IN JUNE 1978 x 
* x 
* MODIFIED TO PERMIT NUMERICAL OPTIMIZATION - 
* USING COPES/CONMIN ® 
* BY LCDR WILLIAM A. DAVIS, JR. IN SEPTEMBER 1980 as 
* MODIFIED TO PERMIT USE OF THE ADS OPTIMIZATION * 
* ROUTINE BY LT. G.L. OLSON IN JUNE 1992 = 
* 
* x 


RAR RRREAKRERKRRKRKRKRRKRRRRRKRKRRRRRRRKRRRRRKRRRRRKRKRRR RRR RRR aR 


CHARACTER*20 NAME 


. COMMON/ADS/DF(21),G(10),IDG(100),IGRAD, INFO, IOPT,IONED,IPRINT, 
sISTRAT, IWK( 2000) ,12(30),0BJ,S(2),VLB(2),VUB(2) ,W( 21,30) ,WK(5000), 
>NCOLA,NCON,NDV,NGT,NRA,NRIWK, NRWK 


COMMON /OLLIE/A( 200,50) ,AMTOT(200),APS,B(3),BFIN,BOA,BVIN,C(3), 
SCANGL ,CF( 200) ~CKpCLI, COF(S), cw(200), DEL(200), DMDOT(200), EA(3, 
-EPS(200), EZERO, F(200, 1),FANGL, FNOBIJ(100), H(200), HINF, HZ(200), 
:QB(200), QINC(200), QTINC(200), OToOT, QTOTAL(100) ,R(200),RB(200), 
: RBASEL, R21, RHOF( 200), ROOTI(4),RPM, ROOTR(4),1T(200), TALFA, TB(200), 
otec(200),TE(200), THICK, THICKI , TIB{ 200) ;TINE,TS(200),TSAT,TSS, 
:TT(200) ,TZ,UF(200),X(200) ,XCOF(5),XPLOT(200),Y(200),2(200),2Z0A, 
sDOBF,DOTH, ICOR(200, 3) MEFE, JINT/JLC,JTC, KFF(50), KFIN(50),KT,NBAN, 
sNEL,NFIN,NSNP 


oe 


GUIDE TO FORTRAN VARIABLE NAMES 


AFOVAS FIN AREA/SMOOTH AREA 


ALFA FIN HALF ANGLE (RADIANS) 

BFIN HEIGHT OF FIN (INCHES) 

BOA TANGENT OF THE FIN HALF ANGLE 
BVIN HEIGHT OF FIN (FEET) 


CALFA COSINE OF ALFA 

CANGL CONE HALF ANGLE (DEGREES) 

CBASE INSIDE CIRCUMFERENCE OF CONDENSER (FEET) 
CExaT INSIDE CIRCUMFERENCE AT CONDENSER EXIT (FEET) 


CE THERMAL CONDUCTIVITY OF CONDENSATE FILM (BTU/HR FT F) 
CL CONDENSER LENGTH (FEET) 

CLI CONDENSER LENGTH (INCHES) 

CPHI COSTNE OF PHT 

CRE CONVERGENCE CRITERION 
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AANAAAAAAANAAAAAAAAAAAQANAAAAAANAAANAAAAAAAQAAAAANAANAAAAQANQAANAANAAANQANAANANAAANAAARAN 


Pcl ea Tet en hn =) 


FILM THICKNESS 

GRADIENT OF OBJECTIVE 

FLOATING POINT VALUE OF NDIV 

CONDENSATE MASS FLOW RATE 

TROUGH WIDTH INCLUDING INCREMENTAL CHANGE 
TROUGH WIDTH AT CONDENSER EXIT 

TROUGH WIDTH AT START OF CONDENSER 

FIN BASE WIDTH 

FORCE VECTOR OF SYSTEM 


"FIN HALF ANGLE (DEGREES) 


CONSTRAINT VALUES ASSOCIATED WITH CURRENT DESIGN 
CONVECTIVE HEAT TRANSFER COEFFICIENT (BTU/HR FT2 F) 
LATENT HEAT OF VAPORIZATION (BTU/LBM) 

CONSTRAINT TYPE IDENTIFIER 

THE ELEMENT NUMBER 

NO. OF ROWS MINUS ONE OF THE UPPER TRIANGULAR FIN 
EQUALS 0 FOR COPPER, AND 1 FOR STAINLESS STEEL 
EQUALS 0 FOR WATER, AND 1 FOR FREON 

GRADIENT CALCULATION IDENTIFIER 

CONTROL PARAMETER 

ONE DIMENSIONAL SEARCH IDENTIFIER 

OPTIMIZER IDENTIFIER 

A FOUR DIGIT PRINT CONTROL 

OPTIMIZATION STRATEGY IDENTIFIER 

NO. OF COLUMNS PLUS ONE BELOW TRIANGULAR FIN 
NUMBER OF SYSTEM NODAL POINT LOCATED AT 

THE CENTER OF SYSTEM COORDINATE 

NUMBER OF SYSTEM NODAL POINT LOCATED AT 

THE JUNCTION OF THE SYMMETRY BOUNDARY AND 

THE LINE OF INTERSECTION BETWEEN THE FIN 

AND THE CONDENSER WALL 

NUMBER OF SYSTEM NODAL POINTS LOCATED ALONG 

THE FIN CONVECTIVE BOUNDARY 

NUMBER OF SYSTEM NODAL POINTS LOCATED ON THE 
SYSMMETRIC BOUNDARY OF TRIANGULAR FIN SECTION 
NOT COUNTING POINTS AT BASE AND APEX 

NUMBER OF COLUMNS WITHIN THE TROUGH WALL SECTION 
SYSTEM BAND WIDTH 

LAST ELEMENT AT BOTTOM SIDE 

FIRST ELEMENT AT BOTTOM SIDE 

NUMBER OF CONSTRAINTS 

NUMBER OF ROWS WITHIN THE FIN 

NUMBER OF ROWS WITHIN THE TROUGH 

NUMBER OF INCREMENT 

NUMBER OF DESIGN VARIABLES 

ELEMENT NUMBER AT BASE OF FIN 

NUMBER OF ELEMENTS 

ELEMENT NUMBER AT END OF TROUGH 

NUMBER OF ROWS IN ARRAY A 

DIMENSIONAL SIZE OF WK 

NUMBER OF SYSTEM NODAL POINTS 

VALUE OF THE OBJECTIVE FUNCTION ASSOCIATED WITH X 
ROTATIONAL SPEED OF HEAT PIPE (RAD/HR) 

CONE HALF ANGLE (RADIANS) 

Pe 

DISTANCE FROM CENTERLINE OF THE HEAT PIPE TO HALF THE F 
HEIGHT 

INSIDE RADIUS OF CONDENSER BASE (FEET) 

INSIDE RADIUS OF CONDENSER BASE (INCHES) 

INSIDE RADIUS OF CONDENSER EXIT (FEET) 
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RPM REVOLUTIONS PER MINUTE 


VECTOR OF DESIGN VARIABLES 


SALFA SINE OF ALFA 

SPHI SINE OF PHI 

SURFAR SURFACE AREA 

rec K CONDENSER WALL THICKNESS (FEET) 
THICKI CONDENSER WALL THICKNESS (INCHES) 


Tel TANGENT OF PHI 
AMBIENT TEMPERATURE 
VISCOSITY 
VLB LOWER BOUNDS ON THE DESIGN VARIABLES ° 
VUB UPPER BOUNDS ON THE DESIGN VARIABLES 
REAL WORK ARRAY 
Z2FIN NUMBER OF FINS 
ZOA RATIO OF TROUGH WIDTH TO FIN BASE WIDTH 


ZSTAR SURFACE LENGTH OF THE FIN MINUS THE SURFACE LENGTH 


COVERED BY THE CONDENSATE IN THE TROUGH 


ZZERO SURFACE LENGTH OF FIN 


PRINT*, ’INPUT FILE NAME’ 
READ*, NAME 

OPEN(10,FILE=NAME) 
OPEN(15,FILE=’/HTPIPE OUTPUT’ ) 
OPEN(14,FILE=’/DUMP OUTPUT’ ) 
OPEN(13,FILE=’/GRAPH OUTPUT’ ) 


THE FOLLOWING READS INPUT DATA, PERFORMS HEAT 
TRANSFER ANALYSIS, AND PRINTS RESULTS. 


ke RR INPUT MODE hake 


ELEMENT CONNECTIVITIES 


READ (10,420) NEL,NSNP,NBAN,IFLUID,IFIN 

WRITE (15,430) NEL,NSNP,NBAN 

Weeree (15,435) ILFLUID,I FIN 

WRITE (15,436) 

WRITE (15,437) 

mene lO,440) {(ICL,(LTCOR(IEL,1),141,3),I1ED*1,NEL) 
WRITE (15,450) 


THE CONDENSER GEOMETRY 
Reo 10,460) CLI,CANGL,RBASEI,R2I,THICKI,BFIN,T2Z 
WRevre, (15,4/70)8 CLI ,CANGL,RBASEL,R2I,THICKI,BFIN,TZ 
READ (10,480) NDIV,NEST,NEFB,NBOTI,NBOTF 
WRITE (15,490) NDIV,NEST,NEFB,NBOTI,NBOTF 

DATA FOR RUNNING 
Reap (10,500) RPM,TSS,TINF,HINF 
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*Q 


*Q 


Q O *+ * 


WRITE (15,510) RPM,TSS,TINF,HINF 
THE CONVERGENCE CRITERIAN 


READ (107520) 9¢F i 
WRITE (157530) 9enEtT 


INTERNAL FIN GEOMETRY 


READ (10,540) FANGL,NFIN 

WRITE (15,550) FANGL,NFIN 

WRITE(*,*) FANGL,NFIN 

READ (10,560) IFF 

WRITE (15,570) IFF 

READ (10,580) (KFIN(I),KFF(I),I#1,IFF) 
READ (10,590) NDOBF,NDOTH,JTC,JLC,JINT,KT 
NHB=NEFB/2 

NBF=NBOTF+1 

DOBF = FLOAT(NDOBF ) 

DOTH = FLOAT(NDOTH) 

WRITE (15,600) ICOR(NBOTI,2),ICOR(NEFB,1),ICOR(NEST,1), 


LICOR(NBOTF,1) 


SET CONSTRAINTS 
NRA=21 

NCOLA=30 

NRWK= 5000 
NRIWK=2000 

NDV= 1 

NCON= 2 

IGRAD= 0 


INITIAL DESIGN 
S(1)= BFIN 


BOUNDS 
VLB(1)= 0.0000001 
VUB(1)= 0.75 


IDENTIFY CONSTRAINTS 
NONLINEAR CONSTRAINT 
IDG(1)=1 

LINEAR CONSTRAINT 
IDG(2)=2 


PRINT*, ‘INPUT THE VALUES FOR ISTRAT,IOPT,IONED AND IPRINT’ 
READ*, ISTRAT,IOPT,  IONED, IPRINT 


INITIALIZE COUNTER 

NO=0 .0 

CHANGE THE INTERNAL PARAMETERS 
INFO=-2 

CALL DADS(INFO,ISTRAT, IOPT, IONED, IPRINT, IGRAD,NDV,NCON,S,VLB, 


>VUB,OBJ,G,IDG,NGT,IZ,DF,W,NRA,NCOLA,WK,NRWK, IWK, NRIWK ) 


IWK(2)=0 
IWK(3)=200 
IWK(5)=4 
IWK(7)=500 
WK(3)=2-0.5 


WK(6)=-0.01 
WK(8)=0.05 
WK(9)=0.50 
WK(10)=0.05 
WK(11)=0.005 
WK(13)=0.001 


WK(14)=0.0001 
WK(21)=0.004 
WK(22)=0.002 
WK(26)=0.004 
WK(37)=0.0000000001 


10 CALL DADS( INFO, ISTRAT,IOPT, IONED, IPRINT, IGRAD,NDV,NCON,S,VLB, 
>VUB,OBJ,G,IDG,NGT,1Z,DF,W,NRA,NCOLA,WK,NRWK, IWK, NRIWK ) 
IF (INFO.EQ.0) GO TO 360 


Rk RRR EXECUTION MODE = ***x** 
NO=NO+1 


CONVERT UNITS OF ALL DIMENSIONAL PARAMETERS 
TO FEET. CONVERT UNITS OF ANGLES TO RADIANS. 
R2I=RBASEI-~-(S(1)/2) 

CL2=CLI/12.0 

R2=R21/12.0 

RBASE=RBASEI/12.0 

BVIN=S(1)/12.0 

DIV=FLOAT(NDIV) 

PI=3.14159265358979 
PHI=2.0*CANGL*PI/360.0 

SPHI=SIN( PHI) 

CPHI=COS( PHI) 

TPHI=TAN( PHI) 

DELX=CL/DIV 

CBASE=2.0*PI*RBASE 
REXIT=RBASE+CL*TPHI 
CEXIT=2.0*PI*REXIT 

THICK=THICKI/12.0 
ALFA=FANGL*2.0*P1I/360.0 
SALFA=SIN(ALFA) 

CALFA=COS(ALFA) 

TALFA=TAN (ALFA) 
EZERO=2.0*(S(1)/12.0)*TALFA 

ZOA=( ( CBASE-( EZERO*NFIN) ) /NFIN) /EZERO 


BOUNDARY CONDITIONS AND TEMPERATURE ESTIMATES 
ALONG THE FIN BOUNDARY 


DO 20 NTINF=NBOTI ,NBOTF 

20 TS(NTINF)=TINF 
DO 30 NNT=NBF,NEL 
TS(NNT)=0.0 

30 H(NNT)=0.0 
DO 40 IGT=1,NEST 
IE=ICOR(IGT,2) 

40 T(IE)=T2Z 
IG=ICOR(NEST,1) 
T(IG)=TZ 
OMEGA IS IN RADIANS/HOUR 
OMEGA=RPM*2.0*PI*60.0 
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QAaAN 


ANaANQ 


50 


60 


65 


70 


80 


DO 50 KL=NBOTI,NBOTF 
H(KL)=HINF 
HIFN=HINF 
TSAT=TSS 
EPSO=ZOA* EZERO 
BOA=TALFA 
SURFAR=NFIN*(2.0*(S(1)/(12*CALFA) )+EPSO) 
EPSEX=(CEXIT-(NFIN*EZERO) ) /NFIN 
BETA=( EPSEX-EPSO) /DIV 
ZZERO=(S(1)/12)/CALFA 
AFOVAS=(ZOA+(1./SALFA))/(1.+Z0A) 
ZA=0.0 
DO 60 NSAT=1,NEST 
TS(NSAT)=TSAT 
TSOLID=(TSAT+TINF) /2.0 
TEMPORARY CHANGE - TFILM 
ASMOOTH=0.0 
ACASE=0.0 
QT=0.0 
OBJ=0.0 
QT1=0.0 
QTF=#0.0 
QTRF=0.0 
QTOT=0.0 
DMTOT=0.0 
NK=NDIV+1 
DO 350 NI=1,NK 
R IS THE INCREMENTAL CHANGE IN THE RADIUS OF THE CONDENSER 
R(NI)=#R2+NI*DELX*SPHI 
RB( NI) =RBASE+NI*DELX*SPHI 
EPS IS THE INCREMENTAL CHANGE IN THE TROUGH WIDTH 
EPS (NI) =*EPSO+NI*BETA 
APS=EPS(NI) 


NODAL POINT COORDINATES 


CALL COORD 

Z(1)=ZA 
DO 70 IZEL#=#=1,NEFB 
NA=ICOR(IZEL,1) 
NB=ICOR(IZEL,2) 
XE=X(NA)-X(NB) 
YE=Y(NA)-Y(NB) 
ELZ=SQRT ( XE**2+YE**2 ) 
Z(IZEL+1)=Z(IZEL)+ELZ 
XZB=X(ICOR(NHB,1))-X(ICOR(1 
YZB=Y(ICOR(NHB,1))-Y(ICOR(1 
ZB=SQRT (XZB**2+YZB** 2 ) 
ZC=ZZERO 

IM=1 


,2)) 
,2)) 


PARABOLIC TEMPERATURE DISTRIBUTION ALONG THE FIN 
BOUNDARY ,USING LAGRANGE INTERPOLATION 


TPlar( LCorn Ge, 2 ))) 
TPZeT(LCOR( NHB, 2 ))) 
TP3=T( ICOR(NEFB,1) ) 
AP1=TP1/(ZB*ZC) 
AP2=TP2/(ZB*( ZB-ZC) ) 
AP3=TP3/(ZC*(ZC-ZB) ) 
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BP1l=-(2ZB+ZC)*AP1 
BP2=-ZC*AP2 
BP3=-ZB*AP3 
A12AP1+AP2+AP3 
B1=BP1+BP2+BP3 
TC=0.0 
DO 90 NY=#1,NEST 
90 TC=TC+T(ICOR(NY,2) ) 
AY=FLOAT(NY+1 ) 
TF=(TC+T(ICOR(NY,1))+AY*TS(NY))/(2.0*AY) 


SOLID-FLUID PROPERTIES 
WATER PROPERTIES 


PEVPPLUID.EQ.1) GO TO 91 

HFG=1093.88-0 .5703*TS(1)+0.00012819%*(TS(1)**2) 
1-0.0000008824*(TS(1)**3) 
RHOF(NI)#62.774-0.00255698*TF-0.000053572*TEF**2 
CF(NI)#0.3034+0.000738927*TF-0 .00000147321*TF**2 
UF(NI)=0.001397-0.000014669*TF+0.0000000631253*TF**2-0.00000 
100000976569*TEX** 3 


FREON PROPERTIES 


91 IF(IFLUID.EQ.0) GO TO 92 
HFG=69 .5459-0.0156011*TS(1)-0.000455294*(TS(1)**2)+0.00000104144* ( 
1TS(1)**3) 
RHOF(NI)#102.059- 0. 025364*TF-0 .000502649*(TF**2)+0.00000135407*(TF 
1**3) 
CF(NI)=0.0594858-0.000429765*TF+0 .00000348218*TF**2-0.000000010416 
L8*TR**3 
UF(NI)=0.00078-0.00000525*TF+0.0000000125*TF**2 
92 UF(NI)=3600*UF(NTI) 
Det PRIN. EO.1) GO TO 93 
CW(NI)#231.7772-0.02222*TSOLID 
93 IF(IFIN.EQ.0) GO TO 94 
CW(NI)=8.776+0. 00265*TSOLID 
94 CK=CW(NI) 
CONST=RHOF (NI) ** 2*OMEGA* *2*HFG*CPHI *CALFA*R(NI) 


INITIAL FILM THICKNESS 
HE {NI .GT.1),GO TO 100 
DEL(1)=21.107%*( ((TSAT-TINF) *CF(NI)/(UF(NI) *HFG) )**. 25)*((UF(NI)/( 
LRHOF (NI) *OMEGA) ) **0.5) 
100 CONTINUE 


AVERAGE ELEMENT CONVECTIVE COEFFICIENT ALONG 
THE FIN BOUNDARY 


ZSTAR=ZZERO-DEL(NI)/CALFA 

AZZ=DEL(NI)/SALFA 

ZZ2=ZSTAR 

HDEN=( (-A1*ZZ**3)/3.0-(B1*ZZ**2)/2.0)+Z22*(TSAT-T(1) ) 
AZS=ABS(4*CF(NI) *UF(NI) *HDEN/CONST) **0.25 

HAC=0.0 

DO 190 IEL=1,NEFB 

AZ=Z(IEL) 

BZ=Z(IEL+1) 

Te (ZSTAR.LE.BZ) GO TO 110 
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AANA 


AAN 


AAD 


GO TO 120 

110 IF (HAC.NE.0.0) GO TO 180 
BZ2ZSTAR 

120 DF (TEL. NE.2) Go TO 130 
AK=(BZ-AZ)/5.0 





ZZ2AK 
GO TO 140 

130 AK=(BZ-AZ)/4.0 
ZZ22AZ 


140 ZEL=4*AK 
DO 150 NH=1,5 
HDEN=(~-1.0*(A1*ZZ**3/3.0+B1*ZZ**2/2.0) )+ZZ*( TSAT-T(1) ) 
HZ(NH)=ABS(CF(NI)**3*CONST/(4*UF(NI) *HDEN) )**0.25 | 

150 22=Z2Z2+AK | 
AVERAGE H USING SIMPSONS RULE 
CONH@AK* (HZ(1)+4*HZ(2)+2*HZ(3)+4*HZ(4)+HZ(5))/(3*ZEL) : 
IF (ZSTAR.EQ.BZ) GO TO 160 
H( IEL)=CONH 
GO TO 190 

160 AZ=ZSTAR 
HAZ=CONH*(AZ—-Z(IEL) ) 

DELA2AZS 
170 BZ=Z(IEL+1) 
DELB=(BZ-ZSTAR) *AZZ/(ZZERO-ZSTAR) 
DELZ=(DELA+DELB) /2.0 
HAC=#(BZ-AZ) *CF(NI) /DELZ 
H( TEL) =(HAZ+HAC) /(BZ-Z(IEL) ) 
GO TO 190 
180 AZ=Z(IEL) 
DELA=DELB 
HAZ=0.0 
GO TO 170 

190 CONTINUE 
NETI=NEFB+1 
DO 200 IEL2NETI,NEST 

200 H(IEL)=CF(NI)/DEL(NI) 


ENTRY INTO THE FINITE ELEMENT SOLUTION 


CALL FORMAF 
CALL BANDEC (NSNP,NBAN,1) 


THE TEMPERATURE DISTRIBUTION 


DO 210 NT=1,NSNP 

210 T(NT)=F(NT,1) 
TIB(NI)=T( ICOR(NBOTI, 2) ) 
TT(NI)=T( ICOR(NEFB,1) ) 
TE(NI)=T(ICOR(NEST,1) ) 
TB(NI)=#T(ICOR(NBOTF,1)) 
TTS=0 .0 
DO 220 NS=1,NSNP 

220 TTS=TTS+T(NS) 
PN#FLOAT(NS) 
TSOLID=TTS/PN 


Q AT THE BOTTOM SIDE 


QBI=0.0 
DO 230 IBEL=NBOTI,NBOTF 
a2 


NKA#ICOR(IBEL,1) 
NKB#=ICOR(IBEL, 2) 
XB=X(NKA) -—X(NKB) 
YB=Y(NKA)-Y(NKB) 
ELB=SQRT( XB**2+YB**2 ) 
230 QBI=QBI+(T(NKA)+T(NKB)-2*TS(IBEL) ) *ELB*H(IBEL)/2.0 
QB(NI)=QBI*DELX 


ITERATION UNTIL CONVERGENCE CRITERIA IS MET 


IF (IM.EQ.1) GO TO 240 
QJ=QBI 
GO TO 250 

240 QI=QBI 
IM=2 
GO TO 80 

250 AQ=#ABS(QJ-QI)/QJ 
IF (AQ.LE.CRIT) GO TO 260 
QI=QJ 
GO TO 80 

260 DMDOT(NI)=2.*QBI*DELX/HFG 
DMTOT=DMTOT+DMDOT( NI) 
Cl=RHOF (NI) **2*OMEGA**2*R(NI) *SPHI/(3*UF(NI) ) 
XCOF(1)=-DMTOT 
XCOF(2)=0.0 
XCOF(3)=0.0 
XCOF(4)=C1*EPS(NTI) 
XCOF(5)=C1*TALFA 
M=4 
CALL DPOLRT (M,IER) 


Pte (ROOTRKS ) GT 20.0) GO TO 270 
TR (ROOTR(2).GT.0.0) GO TO 280 
tee (ROOTR(3).GT.0.0) GO TO 290 
IF (ROOTR(4).GT.0.0) GO TO 300 


THE CONDENSATE THICKNESS 


270 DEL(NI+1)=#ROOTR(1) 
GO TO 310 
280 DEL(NI+1)#ROOTR( 2) 
GO TO 310 
290 DEL(NI+1)=ROOTR( 3) 
GO TO 310 
300 DEL(NI+1)=ROOTR( 4) 
310 QEL=0.0 
TPO Gh .NE.2) sGO.TO 320 


OPFROMe THE TOP SIDE 
Q THROUGH FIN 


QEL=0.0 
320 DO 330 IQEL=1,NEFB 

KA*ICOR(IQEL,1) 

KB=ICOR(IQEL, 2) 

XQEL=X(KA)—-X( KB) 

YQEL=Y(KB)-Y(KA) 

ELM*SQRT( XQEL**2+YQEL**2 ) 

QEL=QEL+(2*TS(IQEL)—-T(KA)-T( KB) ) *ELM*H(IQEL)/2.0 
330 CONTINUE 
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QINC(NI)=QEL*DELX 
AMTOT(NI)=DMTOT 
QET=QEL*DELX*NFIN*2 
QT=QT+QET 

QA=QBI *DELX*NFIN*2 


Q THROUGH TROUGH 


QTRF=0.0 

DO 340 IQEL=NEFB+1,NEST 

KA=ICOR(IQEL,1) 

KB=ICOR(IQEL,2) 

XQEL#=X(KA)=—X( KB) 

YQELY( KB) -Y( KA) 

ELM=SQRT( XQEL**2+YQEL**2) 

QTRF=QTRF+(2*TS(IQEL)-T(KA)—-T( KB) ) *ELM*H(IQEL)/2.0 
340 CONTINUE 

QTINC(NI )=QTRF*DELX 

QTOTAL(NI)=QINC(NI)+QTINC(NI) 

QTRFT=QTRF*DELX*NFIN*2. 

QTF=QTF+QTRFT 

QTOT=QTOT+QA 

ASMOOTH=#2*PI*RB(NI) *DELX+ASMOOTH 

ACASE#ACASE+NFIN*(((2*S(1))/(12*CALFA) )+EPSO) *DELX 
350 CONTINUE 

EVALUATE OBJECTIVE FUNCTIONS AND CONSTRAINTS 

OBJ=—(QTOT) 

WRITE(15,*) ’OBJECTIVE FUNCTION="’, OBJ,’BFIN#=’, S(1) 

FIRST CONSTRAINT IS TO ENSURE THE RATIO ZOA IS NOT NEGATIVE 

G(1)=-( ( (CBASE- (2*NFIN*S(1)*TALPA/12))/NFIN)/(S(1)*TALPA/6) ) 

G(1)=1000.0*G(1) 

THE SECOND CONSTRAINT IS TO ENSURE THE CONDENSATE LEVEL IS NO 

GREATER THAN THE FIN HEIGHT 

G(2)=#=-((8S(1)/12.0)-DEL(NI)) 

XPLOT(NO)#=#S(1) 

FNOBJ (NO) =—-OBJ 

ARATIO#ACASE/ASMOOTH 

WRITE(13,*) XPLOT(NO) , FNOBJ(NO) 

WRITE(14,*) ARATIO, FNOBJ(NO) 

GO TO 10 
360 CONTINUE 


BFIN=S(1) 

kRekeKR* OUTPUT MODE wReRKKK 

WRITE (15,630) 

DO 370 NR=1,NK 
370 WRITE (15,640) NR,QINC(NR) ,QTINC(NR) ,QTOTAL(NR) 

WRITE (15,650) QT,QTF 

WRITE (15,660) CLI,CANGL,RBASEI,R2I,THICKI,BFIN,RPM,TSS,TINF, 

LRIT,FANGL,ZOA,IFF 

WRITE (15,661) AFOVAS 

WRITE (15,670) BOA,ZOA,NFIN,BVIN, SURFAR 

WRITE (15,680) 

DO 380 NP=1,NSNP 

TCC(NP)=.5555555*(T(NP)-32) 
380 WRITE (15,690) NP,X(NP),Y(NP),T(NP),TCC(NP) 

WRITE (15,700) 

DO 390 KKL=1,NBOTF 

NKX=ICOR(KKL,1) 

NKY=ICOR( KKL,2) 
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XP=X(NKX)—-X(NKY) 
YP=Y(NKX)-Y(NKY) 
EXY=SQRT(XP**2+YP**2) 
QEP=*ABS((T(NKX)+T(NKY)-2*TS(KKL) ) *EXY*H(KKL)/2.0) 
QEP=QEP*DELX 
390 WRITE (15,710) KKL,H(KKL),EXY,QEP 
WRITE (05,720) CRIT 
WRITE (15,730) HFG,NFIN,H(NBOTF),TSAT,RPM,QTOT, QT, FANGL 
WRITE(15,734) 
WRITE(15,735) ROOTR(1),ROOTI(1) 
Wem TE(LS,735) ROOTR(2),ROOTI( 2) 
WRITE(15,735,) ROOTR( 3,),, ROOTI (3) 
WRITE(15,735) ROOTR(4),ROOTI(4) 
WRITE (15,740) 
DO 400 NR=1,NK 
400 WRITE (15,750) NR,DEL(NR),QB(NR),AMTOT(NR),TIB(NR),TT(NR),TE(NR), 
1TB(NR) 
WRITE (15,760) 
DO 410 NG=1,NDIV,2 
410 WRITE (15,770) NG,CW(NG),CF(NG) ,RHOF(NG),UF(NG),EPS(NG),R(NG), 
LQINC(NG) 
RETURN 


412 FORMAT (8X,E12.5,8X,E12.5) 

420 FORMAT (515) 

430 FORMAT (/2X,15HNO.OF.ELEMENTS=,15,10X,18HNO.OF.SYSTEM N.P.=, 
115,/,2X,13HNO.OF BANDED=,15) 

435 FORMAT (/2X,’'IFLUID=’,15,10X, 'IFIN=’,1I5) 

436 FORMAT (2X,’IFLUID = 0 FOR WATER, AND 1 FOR FREON’) 

437 FORMAT (2X,’IFIN = 0 FOR COPPER, AND 1 FOR STAINLESS STEEL’) 

440 FORMAT (415) 

450 FORMAT (/2X, 7HELEMENT, 10X, 3HNP1,14X, 3HNP2,15X, 3HNP3) 

460 FORMAT (7G10.5) 

470 FORMAT (4X,5HCLI= ,E12.5,/,4X,7HCANGL= ,E12.5,/,4X,8HRBASEI.=,E12. 
15,/,4X%,5HR2I= ,E£12.5,/,4X,8HTHICKI= ,E12.5,/,4X,6HBFIN= ,E12.5,/,4 
2, 4HTZ=.,E12.5) 

480 FORMAT (515) 

490 FORMAT (4X,6HNDIV= ,110,/,4X,6HNEST= ,110,/,4X,6HNEFB= ,110,/,4X,7 
1LHNBOTI= ,110,/,4X, 7HNBOTF= ,110) 

500 FORMAT (4F10.2) . 

510 FORMAT (4X,5HRPM= ,E£12.5,/,4X,5HTSS= ,E12.5,/,4X,6HTINF= ,E12.5,/, . 
14X,6HHINF= ,E12.5) 

520 FORMAT (G10.9) 

530 FORMAT (4X,6HCRIT= ,E12.5) 

540 FORMAT (2G10.5) 

550 FORMAT (4X,7HFANGL= ,E£12.5,/,4X,6HNFIN= ,I5) 

560 FORMAT (15) | 

570 FORMAT (4X,5HIFF= ,110) 

580 FORMAT (1615) 

590 FORMAT (615) 

600 FORMAT (///5X, 4HTIB=,15,10X, 3HTT=,15,/,5X, 3HTE=,15,/ 
1,6X, 3HTB=,1I5) 

610 FORMAT (//10X,17HCRASH, CRASH, CRASH) 

620 FORMAT (//5X,4(E12.7,3X)) 

630 FORMAT (2X, 7HSTATION, 2X, 4HQFIN,17X, 7HQTROUGH, 15X, 6HQTOTAL ) 

640 FORMAT (4X,15,E12.5,10X,E12.5,10X,E12.5) 

650 FORMAT (//,4X,11HQFIN TOTAL=,E12.5,10X,15HQTROUGH TOTAL= ,£12.5) 

660 FORMAT (/////,4X,5HCLI= ,E12.5,5X,7HCANGL= ,E12.5,/,4X,S8HRBASEI= , 
1£12.5,2X,S5HR2I= ,£12.5,/,4X,8HTHICKI= ,£12.5,2X,6HBFIN= ,E£12.5,/,4 
2X,5HRPM= ,E£12.5,5X,S5HTSS= ,E£12.5,/,4X,6HTINF= ,E12.5,4X,6HHINF= ,E 


oS 


+ + + + 


‘@) 


312.5,/,4X,6HCRIT= ,£12.5,4X, 7HFANGLE= ,E12.5,/,4x%, sHZ0A= 7512 
45HIFF= ,110) 

661 FORMAT(4X,’FIN AREA/SMOOTH AREA=’,E12.5) 

670 FORMAT (1H1,//2X,4HBOA=,G12.5,5X, 4HZOA=,G612.5,5X, 5HNFIN=,15, 
1/, SHBVIN=,G612.5,5X,13HSURFACE AREA=,G12.5) 

680 FORMAT (//5X,2HNP,6X,1HX,12X,1HY,12X,1HT,12X, ZHTC) 

690 FORMAT (/2X,13,3X,4(F10.6,3X) ) 

700 FORMAT (//2X, 2HEL,8X,1HH,11X,9HEL-LENGTH,15X, 4HQ-EL) 

710 FORMAT (/2X%,122, 3x, 50225, 3%, 52.5, Un oe oe) 

720 FORMAT (/2X,22HCONVERGENCE CRITERIAN=,E15.8) 

730 FORMAT (1H ,//,5X,4HHFG=,E12.5,7, 5X, IM BNO.OF (PINS=) none 
16HH-OUT=,E12.5,/,5X, 5HTSAT#,E12.5,/,5X, 4HRPM=,E12.5,/,5X,6HQ 
2E12.5,/,5X,6HQFIN =,E12.5,/,5X, LIHHALF-ANGLE=, F8.3) 

734 FORMAT(/5X,’ROOTS:’,5X,’REAL PARTS’,15X,’IMAGINARY PARTS’ ) 

735 FORMAT(15X,E12.5,15X,E12.5) 

740 FORMAT (1H0,6X,1HJ,4X,14HFILM THICKNESS, 8X, 2HQB,10X,8HMASS-T 
1/, 4X, 38HTIB, OX, 2HTT, 10x, ZHIb, oA, ee 

750 FORMAT (1H ,4X,14,4X,F12.10,4X,F10.4,6X,F9.5,6X,F5.1,/, 

L6XG Fo. Ox E> a leo, bio). ae) 

760 FORMAT (1H0,6X,1HJ,6X,6HK-WALL, 4X,6HK-FILM, 3X, 7HDENSITY,4xX,9 
LFILM,/,6X, 7HEPSILON,5X,6HRADIUS, 5X, 4HQINC) 

770 FORMAT (1H ,4X,14,4X,F/7.3,4X,F6.4,4X,F6.3,4X,F9.7, 
1/;,4K,F9.7, 4K ET 5, oA pe ok, 1 xe) 

END 


THIS SUBROUTINE DEFINES THE POSITIONS OF SYSTEM COORDINATE P 
SUBROUTINE COORD 


COMMON/ADS/DF(21),G(10),1DG(100) , IGRAD, INFO, IOPT, IONED, IPRIN 
> ISTRAT, IWK(2000) ,12(30),0BJ,S(2),VLB(2),VUB(2) ,W(21,30) ,WK(5 
>NCOLA,NCON,NDV,NGT,NRA,NRIWK,NRWK 


COMMON /OLLIE/A( 200,50) ,AMTOT(200),APS,B(3),BFIN,BOA,BVIN,C( 3 
:CANGL,CF(200),CK,CLI,COF(5) ,CW(200) ,DEL(200) ,DMDOT(200) ,EA(3 
:EPS(200),EZERO,F(200,1),FANGL,FNOBJ(100),H(200),HINF,HZ(200) 
:0B(200),QINC(200) ,QTINC( 200) , QTOT, OTOTAL({ 100) ,R( 200) Ret ZzaGe 
>RBASEI,R2I,RHOF(200) ,ROOTI(4),RPM,ROOTR(4),T(200),TALFA,TB(2 
:TCC(200),TE(200), THICK, THICKI, TIB( 200) , TINE, T5(200)) toa nee 
:TT(200) ,TZ,UF(200),X(200) , XCOF(S) ,XPLOT( 200) ,¥(200 ) Zz 20eee 
:DOBF ,DOTH, ICOR(200,3),IFF,JINT,JLC,JTC,KFF(50) ,KFIN(50),KT,N 
>NEL,NFIN,NSNP ; | 


DELH IS THE STANDARD DIVISION OF FIN HEIGHT 
DELH=S(1)/(12*DOBF) 
X(1)=0.0 
Y¥(1)=THICK+(S(1)/12) 
N=1 
DO 20 I=1,IFF 
ICA=KFIN(TI) 
ICB=KFF(TI) 
CBA=FLOAT( ICB-ICA) 
AN=0.0 
DO 10 II=ICA,ICB 
X(II)=X(1)+N*AN*DELH* TALFA/CBA. 
Y(II)=yY(1)-N*DELH 
10 AN=AN+1.0 
20 N=N+1l 
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AN=0.0 

ICD=ICB-ICA+1 

BO S00 aeaJTC,ILC,JINT 

X(J)=X(1) 

Y(J)=(1.0-AN/DOTH) * THICK 

DO 30 JJ#1,ICD 

X(J+JJ)=X(J)+JJ* EZERO/(2*( CBA+1.0) ) 
30 Y(J+JJ)=sY(J) 

JJ=#ICD 

DO 40 K=1,KT 

X( J+JJ+K) =X(J+JI)+K*APS/(2.0*KT) 
40 Y(J+JJ+K)=Y(J) 
SO AN=AN+1.0 

RETURN 

END 


THIS SUBROUTINUE IS USED TO FORMULATE THE FEM EQUATIONS 
SUBROUTINE FORMAF 


COMMON/ADS/DF(21),G(10),IDG(100),IGRAD,INFO,IOPT,IONED,IPRINT, 
one , pw 2000),12(30),0BJ,5(2),VLB(2),VUB(2),W(21,30) ,WK(5000), 
:NCOLA,NCON,NDV,NGT,NRA,NRIWK, NRWK 


COMMON/OLLIE/A( 200,50) ,AMTOT(200),APS,B(3),BFIN,BOA,BVIN,C(3), 
:CANGL,CF(200),CK,CLI,COF(5) ,CW( 200) ,DEL(200) ,DMDOT( 200) ,EA(3,3), 
:EPS(200),EZERO,F(200,1),FANGL, FNOBJ(100),H(200) ,HINF,HZ(200), 
:0B(200),QINC( 200) ,QTINC( 200) ,QTOT, QTOTAL(100),R(200),RB(200), 
:RBASEI,R21I,RHOF(200) ,ROOTI(4),RPM,ROOTR(4),1T(200),TALFA,TB(200), 
eeee coe le2u0), LHeCK, THICKL, TIB( 200) ,; TINE, TS(200) ,TSAT,TSS,, 
ser 200) ,TZ,UFr({ 200) ,X( 200), XCOP(5) ,XPLOT(200),Y¥(200),2(200),Z0A, 
:DOBF, DOTH, ICOR(200,3),IfF,JINT,JLC,JTC,KFF(50),KFIN(50),KT,NBAN, 
>NEL,NFIN,NSNP 


DO 20 N=1,NSNP 
F(N,1)20.0 
DO 10 MA=1,NBAN 

10 A(N,MA)=0.0 

20 CONTINUE 
DO 70 IEL=1,NEL 
IA=ICOR(IEL,1) 
IB=ICOR(IEL,2) 
IC=ICOR(IEL, 3) 
B(1)=Y(IB)-Y(IC) 
B(2)=#Y(IC)-Y(IA) 
B(3)=Y(IA)-Y(IB) 
C(1)=X(IC)-—xX(IB) 
C(2)=#X(IA)-X(IC) 
C(3)=X(IB)—X( IA) 
LENGTH BETWEEN ELEMENT NODES 1 AND 2 
EL=SQRT(C(3) **2+B(3) **2) 
AREA OF TRIANGULAR ELEMENT 
AS=ABS((B(1)*C(2)-B(2)*C(1))/2.0) 
HC=H( IEL)/CK 
BO C0 d=), 3 
JJmICOR(IEL,J) 
DO 50 K=1,3 
KK=ICOR(IEL,K) 
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10 


20 


FORMING THE A MATRIX 
EA(J,K)=(B(J)*B(K)+C(J)*C(K) )/(4*AS ) 
IF (HC.EQ.0.0) GO TO 40 
HEL=HC*EL/6.0 

IF "I. EQ. 3) GOnrOn4e 
TEU(K.EQ.3)@GOnTomae 

IF (J.EQ.K) GO TO 30 
EA(J,K)=EA(J,K)+HEL 

GO TO 40 


EA(J,K)=EA(J,K)+2*HEL 
IFS (CKK.LY.JJ) GOnelogso 
NW=KK-JJ+1l 

A(JJ,NW) =A(JJ,NW)+EA(J,K) 
CONTINUE 

CONTINUE 

FORMING THE F MATRIX 
FE=HC*TS(IEL) *EL/2.0 
F(IA,1)=F(IA,1)+FE 
F(IB,1)=F(1IB,1)+FE 
CONTINUE 

RETURN 

END 


EQUATION SOLVER OF A SYMMETRIC MATRIX THAT HAS BEEN TRANS-— 
FORMED INTO BANDED FORM. 


SUBROUTINE BANDEC (NEQ,MAXB,NVEC) 


COMMON/ADS/DF(21),G(10),IDG(100),IGRAD, INFO, IOPT,IONED,IPRIN 
: ISTRAT, IWK( 2000) ,12(30) ,OBJ,S(2),VLB(2),VUB(2) ,W(21,30),WK(5 
:NCOLA,NCON,NDV,NGT,NRA,NRIWK,NRWK 


COMMON/OLLIE/A( 200,50) ,AMTOT(200),APS,B(3),BFIN,BOA,BVIN,C(3 
:CANGL,CF(200),CK,CLI,COF(5) ,CW( 200) ,DEL( 200) ,DMDOT( 200), EA(3 
:EPS(200),EZERO,F(200,1),FANGL, FNOBJ(100),H(200) ,HINF,HZ(200) 
:QB(200) ,QINC(200) ,QTINC(200) ,QTOT, QTOTAL(100) ,R(200),RB( 200) 
: RBASEI,R21I,RHOF(200) ,ROOTI(4),RPM,ROOTR(4),1T(200),TALFA,TB( 2 
:7TCC(200) , TE( 200), THICK, THICKI,TIB(200) , TINE, 1TS(200)) Tsa27 7a. 
:TT(200),T2,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y{ 200)5 2( 20g 
:DOBF,DOTH, ICOR(200,3),IFF,JINT,JLC, JTC, KFF(50) ,KEIN(50)9 te 
>NEL,NFIN,NSNP 


LOOP=NEQ-1 

DO 20 I=1,LOOP 

MB=I+1 

NB=MINO( I+MAXB-1,NEQ) 

DO 20 J=MB,NB 

L=J+2—-MB 

D=A(I,L)/A(I,1) 

DO 10 MM=1,NVEC 
F(J,MM)=F(J,MM)—-D*F(I,MM) 
MM=MINO ( MAXB—-L+1,NEQ—J+1 ) 
DO 20 K=1,MM 

NN=L+K-1 
A(J,K)=A(J,K)—-D*A(I,NN) 
DO 30 I=1,NVEC 


30 F(NEQ,1I)=F(NEQ,1I)/A(NEQ,1) 


DO 50 I=2,NEQ 
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J=NEQ-I+1 

K=MINO(NEQ-J+1,MAXB) 

DO 50 MM=1,NVEC 

DO 40 L=#2,K 

MB=J+L-l 
F(J,MM)=F(J,MM)-A(J,L)*F(MB,MM) 
F(J,MM)=F(J,MM)/A(J,1) 

RETURN 

END 


SUBROUTINE DPOLRT COMPUTES THE ROOTS OF A REAL 
POLYNOMIAL USING A NEWTON-RAPHSON ITERATIVE 
TECHNIQUE. 


SUBROUTINE DPOLRT (M,IER) 


COMMON/ADS/DF(21),G(10),1IDG(100),IGRAD, INFO, IOPT,IONED,IPRINT, 
: ISTRAT, IWK(2000),12(30),0BJ,S(2),VLB(2),VUB(2) ,W(21,30),WK(5000), 
:NCOLA,NCON,NDV,NGT,NRA,NRIWK, NRWK 


COMMON/OLLIE/A( 200,50) ,AMTOT( 200) ,APS,B(3),BFIN,BOA,BVIN,C(3), 
:CANGL,CF(200),CK,CLI,COF(5),CW(200),DEL(200) ,DMDOT(200),EA(3,3), 
:EPS(200),EZERO,F(200,1),FANGL, FNOBJ(100),H(200),HINF,HZ(200), 
:QB(200) ,QINC(200) ,QTINC( 200) ,QTOT, QTOTAL(100),R(200),RB(200), 
:RBASEI,R21I,RHOF( 200) ,ROOTI(4),RPM,ROOTR(4),T(200),TALFA,TB(200), 
o20¢.200),TE(200), THICK, THICKI,TIB(200) , TINF,TS(200),TSAT,TSS, 
:TT(200),TZ,UF(200),X(200) ,XCOF(5),XPLOT( 200) ,¥(200),2(200),Z0A, 
:DOBF, DOTH, ICOR(200,3),IFF,JINT,JLC,JTC, KFF(50),KFIN(50),KT,NBAN, 
>NEL,NFIN,NSNP 


IFIT=0 

N=M 

IER=0 

IF (XCOF(N+1)) 10,40,10 
mee {N) 20520,60 


SET ERROR CODE TO 1 


IER=1 
RETURN 

SET ERROR CODE TO 4 
IER=4 
GO TO 30 

SET ERROR CODE TO 2 
IER2=2 
GO TO 30 
IF (N-36) 70,70,50 
NX=N 
NXX=N+1 
N2=1 
KJ1=N+1 
DO 80 L=1,KJl 
MT#KJ1-L+1l 
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AaAN AAD -AAND 


aAaAN 


AAA 


aAaAN 


80 


90 


100 


Pee 


PZ 0 


3 0 


140 


150 


160 


170 


COF(MT)=XCOF(L) 
SET INITIAL VALUES 


XO=.00500101 
yo=0.01000101 


ZERO INITIAL VALUE COUNTER 


IN=0 
XX=XO 


INCREMENT INITIAL VALUES AND COUNTER 


XO=—-10.0*YO 
YOx2=-10.0*XX 


SET X AND Y TO CURRENT VALUE 


XX=XO 
YYzYO 
INzIN+1l 
GO TO 120 
IFIT=1l 
XPR=XX 
YPR#YY 


EVALUATE POLYNOMIAL AND DERIVATIVES 


ICT=0 : 
UX=0.0. 

UY=0.0 

V=0.0 

YT=0.0 

XT=1.0 

U=COF(N+1) 

IF (U) 140,270,140 

DO 150 I=1,N 

L=N-I+1 | 
XT2ZeBXX*XT-YY*YT 
YT2Ze@XX* YT+YY*XT 
U=U+COF(L) *XT2 
V=V+COF(L) *YT2 

FI=I 
UX=UX+FI*XT*COF(L) 
UYsUY-FI*YT*COF(L) 
XT=2=XT 2 

YT=YT2 

SUMSQ=UX* UX+UY*UY 

IF (SUMSQ) 160,230,160 
DX=(V*UY-U*UX ) /SUMSQ 
XXeXX+DX 
DY=—-(U*UY+V*UX) /SUMSQ 
YYzYY+DyY 

IF (ABS(DY)+ABS(DX)-1.0E-05) 210,170,170 


STEP ITERATION COUNTER 


ICT=ICT+1 
TF 1ET—500) 307 oe; od 
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220 


230 
240 


250 


260 


270 


280 


290 


300 
310 


B20 


330 


mer (rte) 210,190,210 
IF (IN-5) 100,200,200 


SET ERROR CODE TO 3 


IER=3 

GO TO 30 

DO 220 L#l1,NXX 
MT#KJ1-L+l 
TEMP=XCOF(MT) 
XCOF(MT)=COF(L) 
COF(L)=TEMP 

ITEMP=2N 

N#=NX 

NX#I TEMP 

IF (IFIT) 250,110,250 
IF (IFIT) 240,100,240 
XX#XPR 

YY=YPR 

IFIT=0 


IF (ABS(YY/XX)-1.0E-04) 280,260,260 


ALPHA#XKXX+XX 
SUMSQ@XX* XX+ YY * YY 
N=N-2 

Gorre™ 290 

XX=0 .0 

NX=NX-1 

NXX=NXX-1 

YY=0.0 

SUMSQ=0 .0 
ALPHA#XKX 

N=N-1 
COF(2)=COF(2)+ALPHA*COF(1) 
DO 300 L=#2,N 


COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1) 


ROOTI(N2)#YY 
ROOTR(N2) #XX 

N22N2+1 

IF (SUMSQ) 320,330,320 
YY=x2-YY 

SUMSQ=0.0 

GO TO 310 

IF (N) 30,30,90 

END 
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